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We consider laminar flow of incompressible electrolytes in long, straight channels driven by pres- 
sure and electro-osmosis. We use a Hilbert space eigenfunction expansion to address the general 
problem of an arbitrary cross section and obtain general results in linear-response theory for the mass 
and charge transport coefficients which satisfy Onsager relations. In the limit of non-overlapping 
Debye layers the transport coefficients are simply expressed in terms of parameters of the electrolyte 
as well as the hydraulic radius IZ = 2A/V with A and V being the cross- sectional area and perime- 
ter, respectively. In particular, we consider the limits of thin non- overlapping as well as strongly 
overlapping Debye layers, respectively, and calculate the corrections to the hydrauUc resistance due 
to electro-hydrodynamic interactions. 



I. INTRODUCTION 

Laminar Hagen-Poiseuille and electro-osmotic flows 
are important to microfluidics and a variety of lab-on- 
a-chip applications [J, Q and the rapid development of 
micro and nano fabrication techniques is putting even 
more emphasis on flow in channels with a variety of 
shapes depending on the fabrication technique in use. As 
an example the list of different geometries includes rect- 
angular channels obtained by hot embossing in polymer 
wafers, semi-circular channels in isotropically etched sur- 
faces, triangular channels in KOH-etched silicon crystals, 
Gaussian-shaped channels in laser-ablated polymer films, 
and elliptic channels in stretched soft polymer PDMS de- 
vices [4]. 

In this paper we introduce our recent attempts [H, @] in 
giving a general account for the mass and charge trans- 
port coefficients for an electrolyte in a micro or nanochan- 



nel of arbitrary cross sectional shape. To further moti- 
vate this work we emphasize that the flow of electrolytes 
in the presence of a zeta potential is a scenario of key 
importance to lab-on-a-chip applications involving bio- 




logical liquids/samples in both microfluidic 
nanofluidic channels [10, [HI, E, [H, Q \M, \1 



II. LINEAR-RESPONSE TRANSPORT 
COEFFICIENTS 

The general steady-state flow problem is illustrated in 
Fig. [1] where pressure gradients and electro-osmosis (EG) 
are playing in concert [20]. We consider a long, straight 
channel of length L having a constant cross section Q of 
area A and boundary dft of length V. For many purposes 
it is natural to introduce a single characteristic length 
scale 



piO) =Po + Ap 
V(0) =Vo + AV, 






FIG. 1: A translation invariant channel of arbitrary cross sec- 
tion Q of area A containing an electrolyte driven by a pressure 
gradient —Ap/L and by electro-osmosis through the poten- 
tial gradient —AV/L. The channel wall dO, has the electrical 
potential which induces a thin, charged Debye layer (dark 
gray) that surrounds the charge neutral bulk (light gray). 
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which in the context of hydrodynamics is recognized as 
half the hydraulic diameter. Indeed, for a circle of radius 
R this gives IZ = R. 

The channel contains an incompressible electrolyte, 
which we for simplicity assume to be binary and sym- 
metric, i.e., containing ions of charge +Ze and —Ze and 
equal diffusivities D. The electrolyte has viscosity r^, per- 
mittivity e, Debye screening length X^^ and bulk conduc- 
tivity (Jo = eD/\^j^ and at the boundary dVt it has a zeta 
potential The laminar, steady-state transport of mass 
and charge is driven by a linear pressure drop Ap and 
a linear voltage drop AV. With these definitions flow 
will be in the positive x direction. In the linear-response 
regime the corresponding volume flow rate Q and charge 
current / are related to the driving fields by 



= G 



Ap 
AV 



G 



Gil Gi2 
G21 G22 



(2) 



where, according to Onsager relations [2l|, G is a sym- 
metric, G12 = G21, two-by-two conductance matrix. In 



2 



the following we introduce the characteristic conductance 
elements 



G* = 



^eo ^mig 




(3) 



which is the well-known result for a channel of circular 
cross section of radius R = IZ^ Xd- 

III. SUMMARY OF RESULTS 

In the following we summarize our results for the trans- 
port coefficients accompanied by more heuristic argu- 
ments before we in the subsequent sections offer more de- 
tailed calculations. The upper diagonal element is the hy- 
draulic conductance or inverse hydraulic resistance which 
to a good approximation is given by 



G 



11 



(4) 



While there is no intrinsic length scale influencing Gn, 
the other elements of G depend on the Debye screen- 
ing length Xd- This length can be comparable to 
and even exceed the transverse dimensions in nano- 
channels [TO, [O, [lH ^ which case the off-diagonal ele- 
ments may depend strongly on the actual cross-sectional 
geometry. However, for thin Debye layers with a van- 
ishing overlap the matrix elements G12, G21, and G22 
are independent of the details of the geometry. For a 
free electro-osmotic flow, a constant velocity field v^^ = 
{e(/r])AV/L is established throughout the channel, ex- 
cept for in the thin Debye layer of vanishing width. Hence 

Q 



v^^A and 



G 
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G21 = G*Q, Xd <C 1Z. 



From Ohm's law / = {aoA/L)AV it follows that 



G 



22 



G 



mig5 



Xd <C7^. 



(5a) 



(5b) 



For strongly overlapping Debye layers we shall see that 
in general 



G12 — G; 



21 



8X1 



Xd > 11, 

,2 /\2 



(6a) 



G22 = G^ig + 0(7^7A^)), Xd » 7^. (6b) 

We emphasize that the above results are generally valid 
for symmetric electrolytes as well as for asymmetric elec- 
trolytes. We also note that the expressions agree fully 
with the corresponding limits for a circular cross section 
and the infinite parallel plate system, were explicit solu- 
tions exist in terms of Bessel functions 0, 0] and cosine 
hyperbolic functions [23], respectively. From the corre- 
sponding resistance matrix R = G~^ we get the hydraulic 
resistance 



R 



1 



11 



(7a) 



where (3 = Gi2G2i/(GiiG22) is the Debye-layer correc- 
tion factor to the hydraulic resistance. In the two limits 
we have 



2/-2 



/3: 



8e"C 



X < 



1 



\8XlJ 



Xd<.TI 



Ad > 7^ 



(7b) 



For C, going to zero (3 vanishes and we recover the usual 
result for the hydraulic resistance. 



IV. GOVERNING EQUATIONS 

For the system illustrated in Fig. [H an external pres- 
sure gradient Vp = —{Ap/L)ex and an external electri- 
cal field E = Ebx = {AV/L)ex is applied. There is full 
translation invar iance along the x axis, from which it fol- 
lows that the velocity field is of the form v(r) = v{rj_)ex 
where = yey + ze^. For the equilibrium potential 
and the corresponding charge density we have 0eq(r) = 
(/>eq(i*±) and Peq(^) ~ P^g ('^ ± ) ^ respectively. We will use 
the Dirac bra-ket notation [2J, [25| which is mainly ap- 
preciated by researchers with a background in quantum 
physics, but as we shall see it allows for a very com- 
pact, and in our mind elegant, description of the present 
purely classical transport problem. In the following func- 
tions f{r±) in the domain ft are written as |/) with inner 
products defined by the cross-section integral 



{f\g)= [ dr^f{r^)g{r^). 



(8) 



From the Navier-Stokes equation it follows that the ve- 
locity of the laminar flow is governed by the following 
force balance 







l) + r?V? 



I \ Ay, 



(9) 



where = dy -\- d'^ is the 2D Laplacian and |l) cor- 
responds to the unit function, i.e. ^(r^) = 1. The ffist 
term is the force-density from the pressure gradient, the 
second term is viscous force-density, and the third term 
is force-density transferred to the liquid from the action 
of the electrical field on the electrolyte ions. The equi- 
librium potential \(j)eq) and the charge density |Peq) are 
related by the Poisson equation 



Vi|</>eq) 



--\Pe. 



(10) 



hyd 



The velocity \v) is subject to a no-slip boundary condi- 
tion on dft while the equilibrium potential |(/)eq) equals 
the zeta potential ( on dQ. Obviously, we also need a sta- 
tistical model for the electrolyte, and in the subsequent 
sections we will use the Boltzmann model where the 
equilibrium potential |^eq) is governed by the Poisson- 
Boltzmann equation. However, before turning to a spe- 
cific model we will first derive general results which are 
independent of the description of the electrolyte. 



We first note that because Eq. (|9]) is linear we can 
decompose the velocity as — \vp) ^ I'^eo)^ where \vp^ 
is the Hagen-Poiseuille pressure driven velocity governed 
by 



i> + ^vi|^p>, 



and 



is the electro-osmotic velocity given by 



(C|l)-|0eq)). 



(11) 



(12) 



The latter result is obtained by substituting Eq. (p!Q|) for 
in Eq. The upper diagonal element in G is 

given by Gn = (l\vp) / /S.p which may be parameterized 
according to Eq. (j3|). The upper off-diagonal element 
is given by G12 = (l|'Ueo)/AV and combined with the 
Onsager relation we get 

C?12 = G21 = -^^(1|C - 4>.<,) = -^^(C - '^eq), (13) 

where we have used that = A and introduced the 

average potential (^eq = (</^eq|l)/(l|l)- 

There are two contributions to the lower diagonal el- 
ement G22; one from migration, G^^^ = (l|cr)/L, and 
one from electro-osmotic convection of charge, G22"^ = 



■^eq I 



3)/AV, so that 



G 



22 



^mig 
^22 " 



■^22 



^■(lk)-;^(P:q|C-0eq), (14) 

where the electrical conductivity cr(rx) depends on the 
particular model for the electrolyte. For thin non- 
overlapping Debye layers we note that ^eq — so that 
Eq. ([T3|) reduces to Eq. ([5a|) and, similarly since the in- 
duced charge density is low, Eq. (p!4|) reduces to Eq. (|5bp . 
For strongly overlapping Debye layers the weak screen- 
ing means that ^eq approaches C so that the off-diagonal 
elements G12 = G21 and the G22"^ part of G22 vanish en- 
tirely. In the following we consider a particular model for 
the electrolyte and calculate the asymptotic suppression 
as a function of the Debye screening length Xd- 



V. DEBYE-HUCKEL APPROXIMATION 

Here we will limit ourselves to the Debye-Hiickel ap- 
proximation while more general results beyond that ap- 
proximation can be found in Ref. [6]. In the Debye- 
Hiickel approximation the equilibrium potential |0eq) 
is governed by the linearized Poisson-Boltzmann equa- 
tion [3 



(15) 



where Ad is the Debye screening length which for a sym- 
metric electrolyte is given by 



A 



D 



2{Zefco 



(16) 




FIG. 2: Examples of the first 7 eigenfunctions \il^n) of Eq. ([18} 
with the eigenvalue Hi^ increasing with increasing n. For this 
particular case (kiTV)'^ ^ 5.05 and Ai/A — 0.59 while modes 
with n = 2 and n = 4 will in this case have An = due to 
the symmetry. 



with bulk concentration Cq. The Debye-Hiickel approx- 
imation is valid in the limit Z(e <C ksT where ther- 
mal energy dominates over electrostatic energy. Since we 
consider an open system connected to reservoirs at both 
ends of the channel we are able to define a bulk equi- 
librium concentration in the reservoirs even in the limit 
of strongly overlapping Debye layers inside the channel. 
Thus, strongly overlapping Debye layers do in this case 
not violate the underlying assumptions of the Poisson- 
Boltzmann equation. 



A. Hilbert space formulation 

In order to solve Eqs. ([9]), ([10]), and ([T5|) we wih take 
advantage of the Hilbert space formulation [2^, often 
employed in quantum mechanics [25]. The Hilbert space 
of real functions on Q is defined by the inner product 
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Af/A 


a 


7 


circle 


7? - 5.78"'^ 


- 0.69"'^ 


4:7T 




quarter- circle 


5.08^ 


0.65^ 


29.97^ 


0.93^ 


half- circle 


5.52^ 


0.64^ 


33.17^ 


0.99^ 


ellipse(l:2) 


6.00^ 


0.67^ 


IOtt^ 


1.05^ 


ellipse(l:3) 


6.16^ 


0.62"^ 


407r/3^ 


1.11^ 


ellipse(l:4) 


6.28"^ 


0.58"^ 


177r^ 


1.14^ 


triangle(l:l:l) 


47rV9 ^ 4.3r 
(2^^)^-4.23" 


G/tt^ ^ 0.61^ 


20^3" 5/6 0.83" 


triangle(l:lV2) 


512/97r^ 0.58" 


38.33^ 


0.82^ 


square(l:l) 


7v^/2 4.93" 


64/7r^ ^ 0.66" 


28.45^ 


0.89^ 


rectangle(l:2) 


57rV9 - 5.48" 


64/7r^ :^ 0.66" 


34.98^ 


0.97^ 


rectangle (1:3) 


57rV8 ^ 6.17" 


64/7r^ ^ 0.66" 


45.57"^ 


1.07^ 


rectangle(l:4) 


177rV25 ^ 6.71" 


64/7r^ 0.66" 


56.98^ 


1.14^ 


rectangle(l:oo) 


^7v^ 9.87" 


64/7r^ 0.66" 


oo 


-3/2^ 


pentagon 


5.20^ 


0.67^ 


26.77^ 


0.92^ 


hexagon 


5.36^ 


0.68^ 


26.08^ 


0.94^ 



TABLE I: Central dimensionless parameters for different geometries. "See e.g. [28] for the eigenmodes and eigenspectrum. 
^Here, 71 ^ 2.405 is the first root of the zeroth Bessel function of the first kind. "See e.g. [5] and references therein. ^Data 
obtained by finite-element simulations [2^. "See e.g. [30] for the eigenmodes and eigenspectrum. -^See e.g. [26] for a solution 
of the Poisson equation. 



in Eq. ([8]) and a complete, countable set (IV^n)} of or- 
thonormal basis functions, i.e.. 



(17) 



where Snm is the Kronecker delta. As our basis functions 
we choose the eigenfunctions (IV^n)} of the Helmholtz 
equation with a zero Dirichlet boundary condition on dQ, 



1,2,3,.... 



(18) 



The eigenstates of Eq. (p!8|) are well-known from a variety 
of different physical systems including membrane dynam- 
ics, the acoustics of drums, the single-particle eigenstates 
of 2D quantum dots, and quantized conductance of quan- 
tum wires. Furthermore, with an appropriate re-scaling 
of the Laplacian by IZ or A/V the lowest eigenvalue has 
a modest dependence on the geometry ["sil, |32| . Fig. [2] 
shows as an example the 7 lowest eigenstates l^/^n) in a 
particular geometry. With this complete basis any func- 
tion in the Hilbert space can be written as a linear com- 
bination of basis functions. In the following we write the 
fields as 



n=l 

00 

^eq) = C|l) - ^K\^n), 

n=l 

CO 



(19a) 
(19b) 
(19c) 



The linear problem is now solved by straightforward bra- 
ket manipulations from which we identify the coefficients 



as 



Ap 1 eCAV 



c 



r]L 1^1 r]L 1^ {i^n>^DY 



i + KAd)2 



2 (V^n|l) 



1 + (/^uAd) 



2 • 



(20a) 
(20b) 
(20c) 



B. Transport equations 



The flow rate and the electrical current are conve- 
niently written as 



Q = {i\v), 

I={pl^\v)+(ToE{l\l), 



(21a) 
(21b) 



where the second relation is the linearized Nernst- 
Planck equation with the first term being the convec- 
tion/streaming current while the second is the ohmic cur- 
rent. 



10 



-2 



10^ 



10^ 



FIG. 3: Rescaled off-diagonal transport coefficients versus rescaled Debye- layer thickness in the Debye-Hiickel limit. The solid 
line is the exact result for a circle, Eq. (|27|) . and the dashed line shows Eq. ()6ap . The data points are finite-element simulations 
for different cross sections, see inset. 



C. Transport coefficients 

Substituting Eqs. (|19ap and (|19cp into these expres- 
sions we identify the transport coefficients as 



Gil = G 



hyd 



E 



8 An 



G\2 = 

G21 = 

G22 = 
where 



00 ^ 



An 



(eC)2 A {KnXpf An 



An 



(22a) 
(22b) 
(22c) 
(22d) 

(23) 



is the effective area of the eigenfunction |V^n)- The ra- 
tio An/ A is consequently a measure of the relative area 
occupied by \ilJn) satisfying the sum-rule Xl^i An = A. 
We note that as expected G obeys the Onsager relation 
G12 = G21- Furthermore, using that 



Xd d 



[i^i^nXory 



2 dXDl^if^nXo)^' 



(24) 



we get the following bound between the off-diagonal ele- 
ments G12 = G21 and the lower diagonal element G22, 



G 



22 



^mig 



2\d OXd ' 



(25) 



D. Asymptotics and limiting cases 

1. The geometrical correction factor 

In analogy with Ref. Q we define a geometrical correc- 
tion factor 7 = GJjyj/Gii which from Eq. ()22ap becomes 



_8 An 

^1 {KnTZf A 



{KiTZy A 
8 Ai 



(26) 



Its relation to the dimensionless parameter a in Ref. H 
is 7 = a/{2C) where C — j A is the compactness. As 
we shall see 7 is of the order unity and only weakly de- 
pendent on the geometry so that Eq. (|4]) is a good ap- 
proximation for the general result in Eq. ()22ap . 



2. Non- overlapping, thin Debye layers 

For the off-diagonal elements of G we use that [1 + 
{i^uXd)'^]''^ = 1 + (^[(a^uAd)^]. In Section [Vl] we numer- 
ically justify that the smallest dimensionless eigenvalue 



6 



hzl is of the order 1/7^^, so we may approximate the sum 
by a factor of unity, see Table [H If we furthermore use 
that 7 1 we arrive at Eq. ([5a|) for Xd <^ Tl. These re- 
sults for the off-diagonal elements are fully equivalent to 
the Helmholtz-Smoluchowski result [23]. For G22 we use 
that (A>:nAD)^[l + (/^nAD)^]~^ = OKkuXd)'^], thus we may 
neglect the second term, whereby we arrive at Eq. ([5b|) . 



The transport coefficients in Eqs. ()22ap to (|22dp are 
thus strongly influenced by the ffist eigenmode which 
may be used for approximations and estimates of the 
transport coefficients. As an example the column for 7 is 
well approximated by only including the ffist eigenvalue 
in the summation in Eq. ([26|) . In fact, the approximation 
7 ^ 1 is indeed reasonable. 



3. Strongly overlapping Debye layers 



B. Transport coefficients 



In the case of k^iXd ^ I may use the result [1 + 
(/^uAd)^]""^ = (/^uAd)"^ + 0[{f^n^D)~^] which together 
with J ^ 1 gives Eq. (|6a|) for strongly overlapping Debye 
layers. For G22 we use Eq. (|25]) and arrive at the result 
in Eq. pp]) . 



4. The circular case 



For a circular cross-section it can be shown that [23 



g: 



hjn/Xp) 
hin/Xo) 



(27) 



where In is the nth modified Bessel function of the ffist 
kind, and were we have explicitly introduced the variable 
IZ to emphasize the asymptotic dependence in Eq. (|6a|) 
for strongly overlapping Debye layers. We note that we 
recover the limits in Eqs. ([5a|) and (|6a|) ioi Xd and 
Xd ^ respectively. 



VI. NUMERICAL RESULTS 

A. The Helmholtz basis 

Only few geometries allow analytical solutions of both 
the Helmholtz equation and the Poisson equation. The 
circle is of course among the most well-known solutions 
and the equilateral triangle is another example. However, 
in general the equations have to be solved numerically, 
and for this purpose we have used the commercially avail- 
able finite-element software Comsol Multiphysics [2^. 
Fig. [2] shows the results of finite-element simulations 
for a particular geometry. The ffist eigenstate of the 
Helmholtz equation is in general non-degenerate and 
numbers for a selection of geometries are tabulated in 
Table 1. Note how the different numbers converge when 
going through the regular polygons starting from the 
equilateral triangle through the square, the regular pen- 
tagon, and the regular hexagon to the circle. In gen- 
eral, is of the order 1/7^^, and for relevant high-order 
modes (those with a nonzero An) the eigenvalue is typ- 
ically much larger. Similarly, for the effective area we 
find that Ai/A < 4/7^ 0.69 and consequently we have 
An/ A < 1 - 4/72 - 0.31 for n > 2. 



Our analytical results predict that when going to ei- 
ther of the limits of thin non-overlapping or strongly 
overlapping Debye layers, the transport coefficients to a 
good approximation only depend on the channel geom- 
etry through the hydraulic radius 71. Therefore, when 
plotted against the rescaled Debye length A^/T^, all our 
results should collapse on the same asymptotes in the two 
limits. 

In Fig. [3] we show the results for the off-diagonal co- 
efficients obtained from finite-element simulations in the 
Debye-Hiickel limit for three different channel cross sec- 
tions, namely two parabola shaped channels of aspect 
ratio 1:1 and 1:5, respectively, and a rectangular channel 
of aspect ratio 1:5. In all cases we find excellent agree- 
ment between the numerics and the asymptotic expres- 
sions. For the comparison we have also included exact 
results, Eq. ([27|) . for the circular cross section as well as 
results based on only the ffist eigenvalue in Eq. ()22bp . 
Even though Eq. ([27|) is derived for a circular geome- 
try we find that it also accounts remarkably well for even 
highly non-circular geometries in the intermediate regime 
of weakly overlapping Debye layers. 



VII. CONCLUSION 

We have analyzed the flow of incompressible elec- 
trolytes in long, straight channels driven by pressure 
and electro-osmosis. By using a powerful Hilbert space 
eigenfunction expansion we have been able to address 
the general problem of an arbitrary cross section and 
obtained general results for the hydraulic and electrical 
transport coefficients. Results for strongly overlapping 
and thin, non-overlapping Debye layers are particular 
simple, and from these analytical results we have cal- 
culated the corrections to the hydraulic resistance due to 
electro- hydrodynamic interactions. These analytical re- 
sults reveal that the geometry dependence only appears 
through the hydraulic radius 1Z and the correction factor 
7, as the expressions only depend on the rescaled De- 
bye length Xd/71 and 7 1. Our numerical analysis 
based on finite-element simulations indicates that these 
conclusions are generally valid also for intermediate val- 
ues oi Xd- The present results constitute an important 
step toward circuit analysis [2O, [sslj of complicated mi- 
cro and nanofluidic networks incorporating complicated 
cross-sectional channel geometries. 
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